/*
 * noise_app.cpp
 *
 *  Created on: Jun 11, 2011
 *      Author: nzhao
 */

#include "share/Usrdef.h"
#include "core/globalvar.h"
#include "core/Noise.h"
#include "share/parameters.h"
#include "share/gsl_app.h"
#include "share/MatAlgebra.h"
#include "share/pulseseq.h"

int IS_WARNING = 0;
double MAX_ERROR = 0.0;

gsl_interp_accel *acc;
gsl_spline *spline;
/////////////////////////////////////////////////////////////////////////////
// Calculate phase1

void noise_ensemble::calc_phase1(double * phase)
{
	cout << "calculating phase1..." << endl;
	double * integrate_sin_mat, * integrate_cos_mat;

	integrate_sin_mat = new double [ntime * nfrequency];
	integrate_cos_mat = new double [ntime * nfrequency];
	calc_int_sin_bt_mat(ntime, nfrequency, dt, df, integrate_sin_mat);
	calc_int_cos_bt_mat(ntime, nfrequency, dt, df, integrate_cos_mat);

	calc_quantity(integrate_sin_mat, integrate_cos_mat, phase);
	delete [] integrate_cos_mat;
	delete [] integrate_sin_mat;
}

void noise_ensemble::calc_int_sin_bt_mat(int nt, int nf, double dt, double df, double * int_sin_mat)
{
	int i, j, idx;
	double ti, wj;

	for( i = 0; i < nt; i ++)
	{
		ti = i * dt;
		for(j = 0; j < nf; j ++)
		{
			wj = j * df;
			idx = i * nf+ j;
			int_sin_mat[ idx ] = ( j == 0 ? 0.0 :( 1.0 - harmonic_mat_cos[ idx ] ) / wj);
		}
	}
}

void noise_ensemble::calc_int_cos_bt_mat(int nt, int nf, double dt, double df, double * int_cos_mat)
{
	int i, j, idx;
	double ti, wj;

	for( i = 0; i < nt; i ++)
	{
		ti = i * dt;
		for(j = 0; j < nf; j ++)
		{
			wj = j * df;
			idx = i * nf+ j;
			int_cos_mat[ idx ] = (j == 0 ? ti : harmonic_mat_sin[ idx ] / wj);
		}
	}
}

void noise_ensemble::calc_phase1_gsl_int(double * res, int path_idx)
{
	int i, j, k, npiece;
	double phase_i;
	double x1, x2, dx, epsabs, epsrel, abserr, tval;
	size_t neval;

	cout << "calculate gsl_phase1 ..." << endl;

	gsl_function F;
    F.function = & b_t_val;

    epsabs = 1.0e-3;
    epsrel = 1.0e-5;
    dx = 20.0;

    for(i = 0; i < ntime; i ++)
    {
    	cout << "time " << i << " calculated." << endl;
    	tval = i * dt;
    	npiece = ceil(tval / dx);

    	F.params = &path_idx;

    	res[i] = 0.0;
    	for(k = 0; k < npiece; k ++)
    	{
    		x1 = k*dx;
    		x2 = (k+1)*dx;
    		x2 = (x2 > tval ? tval : x2);

    		gsl_integration_qng (&F, x1, x2,  epsabs,  epsrel, &phase_i, &abserr, &neval);
//    		cout << "function calls : =" << neval << " err = " << abserr << endl;
    		res[i] += phase_i;
    	}
    }
}

double spline_bt(double t, void * param)
{
	return gsl_spline_eval (spline, t, acc);
}
void noise_ensemble::calc_phase1_gsl_int1(double * res, int path_idx)
{
	int i, k, npiece;
	double phase_i;
	double x1, x2, dx, epsabs, epsrel, abserr, tval;
	size_t neval;

//	double * bt_mat_dense = new double [npath * ntime];
//	calc_quantity(harmonic_mat_sin, harmonic_mat_cos, bt_mat_dense);
	generate_interpolation_point(path_idx, 1);

	acc = gsl_interp_accel_alloc ();
	spline = gsl_spline_alloc (gsl_interp_cspline, ninterpolation);

	gsl_spline_init (spline, interpolation_time_list, bt_interpolation, ninterpolation);

	cout << "calculate gsl_phase1_interpolation ..." << endl;

	gsl_function F;
    F.function = & spline_bt;

    epsabs = ABS_EPSLON;//1.0e-4;
    epsrel = 0.0;
    dx = GAUSS_INTERVAL;//1.0;

//    int ntt; double dtt;
//    ntt =100; dtt = MAX_TIME / ntt;
    for(i = 0; i < ntime; i ++)
    {
//    	cout << "time " << i << " calculated." << endl;
    	tval = i * dt;
    	npiece = ceil(tval / dx);

    	res[i] = 0.0;
    	for(k = 0; k < npiece; k ++)
    	{
    		x1 = k*dx;
    		x2 = (k+1)*dx;
    		x2 = (x2 > tval ? tval : x2);

    		gsl_integration_qng (&F, x1, x2,  epsabs,  epsrel, &phase_i, &abserr, &neval);
    		res[i] += phase_i;
    	}
    }

	gsl_spline_free (spline);
	gsl_interp_accel_free (acc);
}
/////////////////////////////////////////////////////////////////////////////
// Calculate phase2
void noise_ensemble::calc_phase2_gsl_int(double * res, int path_idx)
{
	int i, k, npiece;
	double phase_i;
	double x1, x2, dx, epsabs, epsrel, abserr, tval;
	size_t neval;

	generate_interpolation_point(path_idx, 2);

	acc = gsl_interp_accel_alloc ();
	spline = gsl_spline_alloc (gsl_interp_cspline, ninterpolation);

	gsl_spline_init (spline, interpolation_time_list, bt_interpolation, ninterpolation);

	cout << "calculate gsl_phase2_interpolation ..." << endl;

	gsl_function F;
    F.function = & spline_bt;

    epsabs = ABS_EPSLON;//1.0e-4;
    epsrel = 0.0;
    dx = GAUSS_INTERVAL;//1.0;

    for(i = 0; i < ntime; i ++)
    {
//    	cout << "time " << i << " calculated." << endl;
    	tval = i * dt;
    	npiece = ceil(tval / dx);

    	res[i] = 0.0;
    	for(k = 0; k < npiece; k ++)
    	{
    		x1 = k*dx;
    		x2 = (k+1)*dx;
    		x2 = (x2 > tval ? tval : x2);

    		gsl_integration_qng (&F, x1, x2,  epsabs,  epsrel, &phase_i, &abserr, &neval);
    		res[i] += phase_i;
    	}
    }

	gsl_spline_free (spline);
	gsl_interp_accel_free (acc);
}

void noise_ensemble::calc_phase2(double * phase2, int path_idx)
{
	cout << "calculating phase 2..." << endl;
	int j, k, k1;
	double aa_val, bb_val, ab_val;
	double a, b, a1, b1, t, ti, w, b2int, sinwt, coswt, wk, wk1, tval;

	for(j = 0; j < ntime; j++)
	{
		cout << "time = " << j << endl;
		phase2[j] = 0.0;
		for(k1 = 0; k1 < nfrequency; k1++)
		{
			for(k = 0; k < nfrequency; k++)
			{
				calc_aa_bb_val(k1, k, j, &aa_val, & bb_val);
				ab_val = calc_ab(k1, k, j);

//				a = amplitude_mat_alpha[k * npath + path_idx];
//				a1 = amplitude_mat_alpha[k1 * npath + path_idx];
//				b = amplitude_mat_beta[k * npath + path_idx];
//				b1 = amplitude_mat_beta[k1 * npath + path_idx];
				a = amplitude_mat_alpha[path_idx * nfrequency + k];
				a1 = amplitude_mat_alpha[path_idx * nfrequency + k1];
				b = amplitude_mat_beta[path_idx * nfrequency + k];
				b1 = amplitude_mat_beta[path_idx * nfrequency + k1];

				b2int = a*a1*aa_val + b*b1*bb_val + a*b1*ab_val;
				phase2[j] += b2int;
			}
		}
	}
//	char sim_file[500];
//	strcpy(sim_file, WORKSHOP);
//	strcat(sim_file, DEBUG_PATH);
//	strcat(sim_file, "dirPhase_");
//	strcat(sim_file, OUTPUT_FILE);
//
//
//	ofstream sim(sim_file, ios::binary);
//	if(!sim) assert(0);
//
//	sim.write((const char *)&NPATH, sizeof(int));
//	sim.write((const char *)&NTIME, sizeof(int));
//	sim.write((const char *)&MIN_TIME, sizeof(double));
//	sim.write((const char *)&MAX_TIME, sizeof(double));
//
//	sim.write((const char *) phase2, NTIME * NPATH * sizeof(double) );
//
//	sim.close();
}

void noise_ensemble::calc_aa_bb_val(int k1, int k, int ti, double * aa_val, double * bb_val)
{
	double tval, wk, wk1;

	wk = df*k;
	wk1 = df*k1;
	tval = dt*ti;

	if(k1 < k)
	{
		aa_val[0] = 0.0;
		bb_val[0] = 0.0;
	}
	else
	{
		if(k1 == 0 && k ==0)
		{
			aa_val[0] = tval;
			bb_val[0] = 0.0;
		}
		else if(k1 > 0 && k ==0)
		{
			aa_val[0] = 2.0 * harmonic_mat_sin[ti * nfrequency + k1 ] / wk1;
			bb_val[0] = 0.0;
		}
		else if(k == k1 && k != 0)
		{
			aa_val[0] = ( wk*tval + harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k] ) * 0.5 / wk;
			bb_val[0] = ( wk*tval - harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k] ) * 0.5 / wk;
		}
		else
		{
			aa_val[0] = 2.0* (wk * harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k1]
			            -wk1 * harmonic_mat_sin[ti * nfrequency + k1] * harmonic_mat_cos[ti * nfrequency + k])
			           / (wk*wk - wk1*wk1);

			bb_val[0] = 2.0* (wk1 * harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k1]
			           -wk * harmonic_mat_sin[ti * nfrequency + k1] * harmonic_mat_cos[ti * nfrequency + k])
			           / (wk*wk - wk1*wk1);
		}
	}
}
double noise_ensemble::calc_ab(int k1, int k, int ti)
{
	double res, tval, wk, wk1;

	wk = df*k;
	wk1 = df*k1;
	tval = dt*ti;

	if(k1 == 0 )
	{
		res = 0.0;
	}
	else if(k == k1 && k != 0)
	{
		res =  harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_sin[ti * nfrequency + k] / wk;
	}
	else
	{
		res = 2.0* (wk * harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_sin[ti * nfrequency + k1]
		      +wk1 * harmonic_mat_cos[ti * nfrequency + k1] * harmonic_mat_cos[ti * nfrequency + k] - wk1)
		      / (wk*wk - wk1*wk1);
	}
	return res;
}


void noise_ensemble::calc_phase2_old(double * phase2)
{
	cout << "calculating phase 2..." << endl;

	int i,j,k, dim_freq, dim_freq2;
	dim_freq = (nfrequency + 1)  * nfrequency / 2;
	dim_freq2= nfrequency * nfrequency;

	gen_bilinear_mat();

	cout << "Bilinear mat generated. " << endl;

	double * mat_aa, * mat_bb, * mat_ab;
	mat_aa = new double [dim_freq * ntime];
	mat_bb = new double [dim_freq * ntime];
	mat_ab = new double [dim_freq2 * ntime];

	cout << "mat_aa/bb/ab generated. " << endl;

	for(i = 0; i < ntime; i ++)
	{
		calc_aa_bb_mat(mat_aa + i * dim_freq, mat_bb + i * dim_freq, i);
		calc_ab_mat(mat_ab + i * dim_freq2, i);
	}
//	char amp_file[200];
//
//	strcpy(amp_file, WORKSHOP);
//	strcat(amp_file, DEBUG_PATH);
//	strcat(amp_file, "mat_");
//	strcat(amp_file, PARA_CHAR);
//	strcat(amp_file, ".dat");
//
//	ofstream amp(amp_file);
//	if(!amp) assert(0);
//
//
//	amp << setw(12) << setiosflags(ios::scientific) << setprecision(10);
//
//	for(i = 0; i < dim_freq* ntime; i++)
//	{
//		amp << mat_bb[i]  << endl;
//		cout<< mat_ab[i]  << endl;
//	}


//	cout << "bb" << endl;
//	for(i = 0; i < dim_freq* ntime; i++)
//	{
//		cout << mat_aa[i] << "," << endl;
//	}
//
//	for(i = 0; i < dim_freq2* ntime; i++)
//	{
//		cout << i <<":ab="  << mat_ab[i] << endl;
//	}

//	amp.close();

	MKL_real_mat_mult(npath, ntime, dim_freq, NO_TRANS, NO_TRANS, CONST_ONE, bilinear_alpha, mat_aa, CONST_ZERO, phase2);
	MKL_real_mat_mult(npath, ntime, dim_freq, NO_TRANS, NO_TRANS, CONST_ONE, bilinear_beta, mat_bb, CONST_ONE, phase2);
	MKL_real_mat_mult(npath, ntime, dim_freq2, NO_TRANS, NO_TRANS, CONST_ONE, bilinear_alpha_beta, mat_ab, CONST_ONE, phase2);

	delete [] mat_aa;
	delete [] mat_bb;
	delete [] mat_ab;

	char sim_file[500];
	strcpy(sim_file, WORKSHOP);
	strcat(sim_file, DEBUG_PATH);
	strcat(sim_file, "dirPhase_");
	strcat(sim_file, OUTPUT_FILE);


	ofstream sim(sim_file, ios::binary);
	if(!sim) assert(0);

	sim.write((const char *)&NPATH, sizeof(int));
	sim.write((const char *)&NTIME, sizeof(int));
	sim.write((const char *)&MIN_TIME, sizeof(double));
	sim.write((const char *)&MAX_TIME, sizeof(double));

	sim.write((const char *) phase2, NTIME * NPATH * sizeof(double) );

	sim.close();
}

void noise_ensemble::gen_bilinear_mat()
{
	int i, j, n, idx, idxfreq, idx_i, idx_j, dim_freq, dim_freq2;
	dim_freq = (nfrequency + 1)  * nfrequency / 2;
	dim_freq2= nfrequency * nfrequency;

	bilinear_alpha = new double [npath * dim_freq];
	bilinear_beta= new double [npath * dim_freq];
	bilinear_alpha_beta = new double [npath * dim_freq2];


	for(n = 0; n < npath; n++)
	{
		idxfreq = 0;
		for(i = 0; i < nfrequency; i++)
		{
			for(j = 0; j <= i; j ++)
			{
				idx = idxfreq * npath + n;
				idx_i = i * npath + n;
				idx_j = j * npath + n;

				bilinear_alpha[idx] = amplitude_mat_alpha[idx_i] * amplitude_mat_alpha[idx_j];
				bilinear_beta[idx] = amplitude_mat_beta[idx_i] * amplitude_mat_beta[idx_j];
				idxfreq ++;
			}
		}
	}

	for(n = 0; n < npath; n++)
	{
		idxfreq = 0;
		for(i = 0; i < nfrequency; i++)
		{
			for(j = 0; j < nfrequency; j ++)
			{
				idx = idxfreq * npath + n;
				idx_i = i * npath + n;
				idx_j = j * npath + n;

				bilinear_alpha_beta[idx] = amplitude_mat_beta[idx_i] * amplitude_mat_alpha[idx_j];
				idxfreq ++;
			}
		}
	}
}

void noise_ensemble::calc_aa_bb_mat(double * aa_mat, double * bb_mat, int ti)
{
	int k1, k, idx;
	double wk1, wk, tval;


		tval = ti * dt;

		idx = 0;
		for(k1 = 0; k1 < nfrequency; k1++)
		{
			wk1 = k1 * df;
			for(k = 0; k <= k1; k++)
			{
				wk = k * df;

				if(k1 == 0 && k ==0)
				{
					aa_mat[idx] = tval;
					bb_mat[idx] = 0.0;
				}
				else if(k1 > 0 && k ==0)
				{
					aa_mat[idx] = 2.0 * harmonic_mat_sin[ti * nfrequency + k1 ] / wk1;
					bb_mat[idx] = 0.0;
				}
				else if(k == k1 && k != 0)
				{
					aa_mat[idx] = ( wk*tval + harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k] ) * 0.5 / wk;
					bb_mat[idx] = ( wk*tval - harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k] ) * 0.5 / wk;
				}
				else
				{
					aa_mat[idx] = 2.0* (wk * harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k1]
					             -wk1 * harmonic_mat_sin[ti * nfrequency + k1] * harmonic_mat_cos[ti * nfrequency + k])
					             / (wk*wk - wk1*wk1);

					bb_mat[idx] = 2.0* (wk1 * harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_cos[ti * nfrequency + k1]
					              -wk * harmonic_mat_sin[ti * nfrequency + k1] * harmonic_mat_cos[ti * nfrequency + k])
										             / (wk*wk - wk1*wk1);
				}
				idx ++;
			}
		}

}

void noise_ensemble::calc_ab_mat(double * ab_mat, int ti)
{
	int k1, k, idx;
	double wk1, wk, tval;

	tval = ti * dt;

	idx = 0;
	for(k1 = 0; k1 < nfrequency; k1++)
	{
		wk1 = k1 * df;
		for(k = 0; k < nfrequency; k++)
		{
			wk = k * df;

			if(k1 == 0 )
			{
				ab_mat[idx] = 0.0;
			}
			else if(k == k1 && k != 0)
			{
				ab_mat[idx] =  harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_sin[ti * nfrequency + k] / wk;
			}
			else
			{
				ab_mat[idx] = 2.0* (wk * harmonic_mat_sin[ti * nfrequency + k] * harmonic_mat_sin[ti * nfrequency + k1]
				      +wk1 * harmonic_mat_cos[ti * nfrequency + k1] * harmonic_mat_cos[ti * nfrequency + k] - wk1)
				      / (wk*wk - wk1*wk1);
			}
			idx ++;
		}
	}
}

/////////////////////////////////////////////////////////////////////
// full energy fluctuation
void noise_ensemble::calc_phaseFull_gsl_int(double * res, int path_idx)
{
	int i, k, q, npiece;
	double phase_i;
	double x1, x2, dx, epsabs, epsrel, abserr, tval, t_q, t_q_1, res_q, sign;
	size_t neval;
	int status;

	generate_interpolation_point(path_idx, 3);

	acc = gsl_interp_accel_alloc ();
	spline = gsl_spline_alloc (gsl_interp_cspline, ninterpolation);

	gsl_spline_init (spline, interpolation_time_list, bt_interpolation, ninterpolation);

	gsl_function F;
    F.function = & spline_bt;

    epsabs = ABS_EPSLON;//1.0e-4;
    epsrel = 0.0;
    dx = GAUSS_INTERVAL;//1.0;

    gsl_set_error_handler_off();
    for(i = 0; i < ntime; i ++)
    {
    	tval = i * dt;
    	res[i] = 0.0;
    	sign = -1.0;

    	for(q = 1; q <= PULSE_NUM + 1; q ++)
    	{
    		sign *= -1.0;
    		t_q = (q == PULSE_NUM + 1 ? tval : tval * getPulse(q));
    		t_q_1 = (q == 1 ? 0.0 : tval * getPulse(q - 1) );
    		res_q = 0.0;

    		npiece = ceil( (t_q - t_q_1) / dx);
    		for(k = 0; k < npiece; k ++)
    		{
    			x1 = t_q_1 + k*dx;
    			x2 = (x1 + dx < t_q ? x1 + dx : t_q);

    			status = gsl_integration_qng (&F, x1, x2,  epsabs,  epsrel, &phase_i, &abserr, &neval);
    			if (status == GSL_ETOL )
    			{
    				IS_WARNING ++;
    				if (fabs(abserr) > MAX_ERROR)
    					MAX_ERROR = fabs(abserr);
    				cout << "Warning::" << abserr << " exceed the threashold ::"
    						<< " path " << path_idx
    						<< " tval " << tval
    						<< " pulse "<< t_q_1
    						<< " picece " << k << endl;
    			}
    			res_q += phase_i;
    		}
    		res[i] += sign * res_q;
    	}
    }

	gsl_spline_free (spline);
	gsl_interp_accel_free (acc);
}
